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We discuss the reflection of light by a rheoscopic fluid (a suspension of microscopic 
rod-like crystals) in a steady two-dimensional flow. This is determined by an order pa- 
rameter which is a non-oriented vector, obtained by averaging solutions of a nonlinear 
equation containing the strain rate of the fluid flow. Exact solutions of this equation are 
obtained from solutions of a linear equation which are analogous to Bloch bands for a 
one-dimensional Schrodinger equation with a periodic potential. On some contours of 
the stream function, the order parameter approaches a limit, and on others it depends 
increasingly sensitively upon position. However, in the long-time limit a local average of 
the order parameter is a smooth function of position in both cases. We analyse the topol- 
ogy of the order parameter and the structure of the generic zeros of the order parameter 
field. 



1. Introduction 

Rheoscopic fluids are suspensions of microscopic rod-like crystals which are brought 
into alignment by a fluid flow. They enable the fl ow to be visualised due to the angular 
dependence of light reflection from the cryst als (Matisse fc Gormanlll984 ) and similar 



suspensions are used to make art installations ( Reedll2002[) and to enhance the appearance 
of cosmetic products (Preston 1984). The patterns produced by rheoscopic agents bear 
a complex relation to the underlying flow, which is not yet thoroughly understood. For 
example, a simple stirring motion produces an increasing tightly wound spiral pattern, 
illustrated in figure [1] (which shows additive mixing of light scattered from red, green 
and blue sources, as illustrated in figure [2]). A priori, it is not clear how the orientation 
of the rods should depend upon position within this pattern, or how it will evolve in 
the long-time limit. In this paper we motivate the definition of an order parameter for 
the alignment of the crystals in two-dimensional flows, and show how it may be related 
to the colour of the reflected light. We use this to analyse the long-time limit of the 
orientation patterns formed by steady flows in two dimensions, where trajectories of the 
small crystals follow contours of the stream function, ip(x,y). 

The motion of sm a ll ellip soidal bodies in steady two-dimensional flows was previously 
considered bv ISzeri ( 1993f) . who showed that there may exist orbits (closed contours 



of il}{x,y)) where the principal axis approaches a constant direction, as well as orbits 
where the axis tumbles and where spiral patterns such as that in figure [T] are seen. He 
also showed that the regions where alignment occurs are characterised by a topologi- 
cal index, which was termed the 'flip number', but which is in fact equal to twice the 
Poincare index (the definition of the Poincare index is illustrated in figure [3]). His results 
raise a variety of interesting questions concerning the textures of rheoscopic flows, which 
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Figure 1. Illustrates a spiral pattern which may be generated by motion of a rheoscopic fluid 
in a two-dimensional cellular flow, with stream function '4)(x,y) — sin(a;) sin(i/)/27r. The flow 
is visualised by reflected light from three different coloured sources, as shown schematically in 
flgure[2]and as described in detail in section |3l The arms of the spiral tighten as time increases. 
The times are (successively from left to right) t = 5, t = 9, and t — 18, and the parameters in 
(pT]) . ([2:21 are ai = 0.95, Q2 = 0.05. 




Figure 2. a The direction and degree of ordering of the axes of the crystals in a rheoscopic 
fluid can be revealed by scattering light from red, green and blue sources arranged around the 
sample, b The degree of order of the particles is described by an order parameter vector ^ lying 
within a unit circle, which points in the predominant direction of alignment, with the magnitude 
^ Id ^ 1 indicating the degree of alignment, c The colour of the scattered light is a function 
of the order parameter: because the orientation of the vector <^ is irrelevant, this colour map is 
symmetric under reflection. 

are resolved in this paper. What is the actual appearance of the system under reflected 
light? In the regions where the crystals tumble, they may still have a preferred alignment. 
Does their alignment approach a time-independent limit in the regions where the crystals 
tumble, and if so, how is this limit approached? Is there an abrupt change in appear- 
ance on crossing from the tumbling region to the aligning region? These questions are 
most directly addressed by analysing an order parameter vector (^{r,t) which describes 
the predominant direction of alignment of the crystals, even in regions where they are 
tumbling. 

We show that in the tumbling regions the order parameter forms a progressively more 
tightly-wound spiral pattern, having an increasingly sensitive dependence upon position. 
In the long-time limit the order parameter will fluctuate on a scale which is below the 
resolving power of the eye, and it is necessary to consider a local average of the ori- 
entations within a small disc, from which we determine an averaged order parameter, 
denoted by (C) . We show how the smoothly- varying local average is calculated, so that in 
the long-time limit the ordering of the rods is described by a smoothly varying function 
(Oil"), defined in both the tumbling and the aligning regions. We find that this function 
has no discontinuity at the boundary between the regions. This raises a further question. 
It is observed that the aligning regions have different Poincare indices, which implies that 
(^)(r) must have some form of singularity in the tumbling region. What is the form of 




Figure 3. The Poincare index is a topological invariant. For a vector field in the plane, the 
Poincare index of a closed curve is the number of 2n clockwise rotations of the vector field as 
the curve is traversed, also clockwise. Curves with a non-zero Poincare index must encircle a 
singularity of the field. Because the axial vector of the rod-like crystals is non-oriented, singular- 
ities with half-integer Poincare index are possible: a is a vortex, b is a core and c is a delta, with 
indices -1-1, i, — i respectively. If the curve encloses two singularities, their indices are added: 
for example the Poincare index of the curve in d is ^ — | = 0. 



these singularities? We identify the normal forms for generic singularities which are nodal 
points of the field {0- These singularities have Poincare index ±i and have structures 
which are analogous to singularities which are seen in the ridge patterns of fingerprints. 

The local average order parameter (^) in the long-time limit is illustrated in figure SI 
for a 'journal bearing' flow (that is, a two-dimensional flow between two non-slip non- 
concentric rotating boundaries) , using the visualisation method illustrated in figure [21 
The details of this example will be discussed in section [H but this figure indicates that 
in some cases the solu tion of this problem may be very complicated. 

In two recent works (jWil kinson. Bezuglvv &: Mehlid2008t Bezuglyy. Wilkinson fcMehlig 
20091) . we have considered the alignment of rheoscopic fluids in response to generic, time- 



dependent flows. If a time-dependent flow does not decay (for example, if the fluid is 
continuously stirred), then there is a usually a positive Lyapunov exponent (meaning 
that infinitesimal separations between fluid elements grow exponentially). In the case of 
random flows we also find singularities which are related to those in fingerprints. How- 
ever, there is an important distinction. In the case we consider here the fingerprint-like 
singularities only emerge after performing a local average of the order parameter in the 
long-time limit. In a random fiow, by contrast, singularities may be observed at short 
times. The reasons for the difference are explained at the end of this paper. 

We assume for simplicity that the crystals are axisymmctric. Our discussion applies 
in the case where the bodies are much smaller than any length scale of the flow, and 
for this reason we describe them as 'particles' for the remainder of this paper. For the 
case of a steady two-dimensional flow, the limit of inflnite aspect ratio is a singular 
case, and we retain the aspect ratio (3 of the particles as a parameter in our equations. 
The equation of motion for the unit vector n aligned with the symmetry axis of the 
microscopic axisymmctric particles was originally obtained by lieffery {19^i_ An ele- 
gant solution of this equation of motion was subsequently given by ISzeri ( 1993f ). who 
showed how the solution of this non-linear equation may be obtained by normalising 
a vector which evolves according to a companion linear equation (Szeri gives credit to 
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Figure 4. Illustrating the reflection of red, green and blue light from rheoscopic fluid in a 
'journal bearing' flow in the long-time limit. This is a steady two-dimensional flow between two 
non-concentric rotating non-slip circular walls. In this example both walls rotate in the same 
direction, with the angular velocity of the inner boundary exceeding that of the outer boundary 
by a factor of 20. The aspect ratio of the elliptical particles is /3 = vTQ. 



earlier works by iBretherton ( 19621 ) and Lipscomb et a\ ( 1988 ). but tliese do not contain 
the general solution). Most of the other literature has applied concepts from dynamical 
syst ems theory to the non-linear system of equations obtained by Jeffery (see, for exam- 
ple ("Shin fc Maxev"l99l':'Szeri. Wiggins fc Leal"l99l':'MaUier fc Maxevlll99l HSzerifcLea] 
[1994; Shin & Maxcy 1997; Ga uthicr, Go ndorct & Rabaud 1998|)), and Szeri's own paper 
makes very limited use of his general solution. 

In this work we combine Szeri's solution with the insight that comes from the analogy 
between the companion linear equation and the time-independent Schrodinger equation 
in one dimension. Typically, the contours of the stream function are closed curves, so 
that the trajectory of a fluid element is periodic in time, with a period T which depends 
upon the contour. The evolution of the companion linear equation for a trajectory on 
a closed contour is analogous to the propagation of t he solution o f a time-independent 
Schrodinger equation in a spatially periodic potential (jZiman Il976l) . The solution of this 
latter problem has bands of energy (the Block bands) where generalised eigenstates exist 
(which take the form of Block waves) , interspersed by band gaps, intervals of energy for 
which the electron cannot propagate. The regions where the particles approach a constant 
direction correspond to the band gaps in the solution of the Schrodinger equation, and 
the Bloch bands correspond to the regions where the particles tumble. In the following we 
refer (for reasons which will be discussed in section [5|) to the regions where particles align 
as kyperbolic bands, and the regions where they tumble as elliptic bands. The analogy 
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with sohd-state physics wiU be useful to readers who are familiar with that field, but it 
is not essential to understanding the paper. 

In section [2] we discuss the equation of motion for the direction vector of a single parti- 
cle. We describe its general solution, and also consider the instructive special case of flows 
with a uniform velocity gradient. In section [3] we consider the order parameter, explain- 
ing th e motivation for the definition which was used by iBezuglvv. Wilkinson fcMehlie 
(|2009l) . and explaining how the order parameter is related to the reflection of light from 
the system. Section 2] discusses the case of recirculating flows in two dimensions, in- 
cluding an a nalogy with B l och's theorem a nd id eas related to Anderson localisation 
(described bv lZimanI (|l976l ): iMott fc Twosd (|l96ll )V We show that the order parameter 
has a quasiperiodic structure in regions where the rods tumble. Section [5] considers the 
calculation of the locally averaged order parameter in the long-time limit, and consid- 
ers topological aspects of the solution. We show that the Poincare index can change by 
at most ±^ on crossing a band (with simple annular topology) where the rods tumble, 
and give a criterion for determining the change of Poincare index. If the Poincare index 
changes upon crossing an elliptic band, there must be a singularity of the average or- 
der parameter. In section [5] we identify normal forms for these singularities. Section [6] 
contains some concluding remarks, and a discuss ion of how the results of this paper dif - 
fer from the case of random fl o ws, co nsidered by IWilkinson. Bezuglvv fc Mehlig ( 2008h ; 
Bezuglvv. Wilkinson fcMehlid (|2009f ). 



2. Equation of motion and its solution 

We consider the motion of very small rigid axisymmetric particles immersed in a fluid 
flow with velocity field v{r,t). We assume they are small compared to the characteristic 
length scale of the flow, and also sufficiently small that they do not interact with each 
other or perturb the velocity field. They are assumed to have negligible inertia so that 
their motion is dominated by viscous forces. The equation of motion of the centre of the 
particle r{t) is then the advective equation r = v{r,t) (we shall use dots to indicate time 
derivatives). For sufficiently small particles, the equation of motion for n{t) can only 
depend upon the gradient of the velocity field, described by a tensor A(i) with matrix 
elements Aij{t) = dvi/drj{r{t),t), where r{t) is the particle trajectory. The motion of 
a unit vector n aligned with the axis of sy mmetry i s dete rmined by the condition that 
the torque on the particle is equal to zero. Ijefferv ( 1922) determin ed the equation of 



motion for n{t) in the case of a spheroidal particle, and Brethertonl (1962.) showed that 
the equation of motion for a general axisymmetric particle is of the same form. The 
equation of motion may be written as 

— =Bn-n(n-Bn) (2.1) 

dt ^ ^ ^ ' 

where B is a matrix derived from the rate of strain matrix A as follows 

B = aiA-a2A'^, ai+a2 = l. (2.2) 

Here A""" is the transpose of the matrix A, and ai, a2 are dimensionless parameters which 
are determined by the aspect ratio of the particle: for an ellipsoid of aspect ratio (3 (with 
/3 ^ 1), Jeffery showed that ai = /3^/(/3^ + l)i ct2 = l/(/3^ + !)■ Jeffery originally wrote 
his equations of motion i n component form, how ever our (j2.ip . (|2.2p are equivalent to 
e quations (10) a nd (12) in lMallier fc MaxevI (|l99ll ) with E = i(A-FA'^) andw = VAtt. 
iJefferv (|l922h also discussed the particular case of the motion of an ellipsoid of revo- 



lution in a uniform shear flow, and showed that (except for the limiting case of a rod) 
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the p a rticle exhibits a tumbUng motion, which has been observed experimentaUy by 
Savad (1985). This has been extended to consider the tumbhng motion of particles in 



cel lular flows ([Mallier fc Maxev 



19911 ). The case of a general linear flow was discussed 



by ISzeri. W iggins fc Leal ( 1991 ). who also discussed the response to a more general 

flow field in the lang u age of dynamical s ystems theory. Most o ther works (for example, 

IShin fc Ma^(|l991f ): [SzerifcLeall(|l994 ): IShin fc Maxevl|l997l ): lGauthier. Gondoret fc Rabaud 
(|l998h ) have also used a dynamical systems approach based upon the nonlinear equation 
of motion for n. 

We can, however, solve (j2.ip in terms of the solution of an auxiliary problem, which is 
linear. Specifically, we solve the equation 



dd 

'dt 



= B{t)d 



(2.3) 



to determine a vector d{t). Here B(t) is the matrix B evaluated at the position reached 
by the particle at time t, that is B(t) — aiA{r{t),t) — a2A'^(r(t), i). This equation is 
solved with the initial condition <i(to) = no, where ng is the initial orientation of the 
particle at time to. Now multiply d{t) by a scalar //(t), chosen such that n{t) — ii{t)d{t) 
is a unit vector. We find that this normalised vector does indeed satisfy equation p.ip . 
We therefore have a solution of the nonlinear equation for n{t) in the form 



n{t) = 



d{t) 
\dit)\ 



(2.4) 



This is an exact and completely general solution for the orientation, in terms of the solu- 
tion of a companion linear problem, equation ()2.3p . Because of the superposition principle, 
it is almost always much easier to analyse a linear problem, even in circumstances where 
exact solutions are not available. We ex ploit this adv antage in the remainder of this pa- 
per. This solution was first obtained bv ISzeri ( 1993 ). but remarkably most subsequent 
papers did not make use of this powerful result. 

Solving equation (j2.4p is sufficient for determining the motion of a single particle with 
a specified initial orientation, but in many cases we wish to consider the motion of many 
small particles, or to obtain the solution for an arbitrary initial orientation. In this more 
general context, instead of solving (|2.4p we determine a matrix M(i) which is the solution 
of 

= B(r(t),t)M (2.5) 

with initial condition M(0) — I (the identity matrix). Given this matrix, the solution 
of (|2.5p is d{t) = M(<)do, for any choice of do, so that a single solution suffices for 
all initial directions. A further generalisation is to consider an arbitrary initial position 
for the particle at time to. Let M(r,i,to) be the solution of ()2.ip for a particle which 
reaches position r at time t, having started at rg at time to. In this most general case 
the orientation is a vector field, n(r,t), and our exact solution becomes: 



n{r,t) 



M{r,t,to)n{ro,to) 
\M(r,t,to)n{ro,to)\ 



(2.6) 



In the case of rod-like particles, where /? oo, the matrix B(i) is equal to the velocity- 
gradient matrix A{t), with elements Aij = dvi/drj. In this case the matrix M(f) has a 
simple physical interpretation, and in the following we use ]V[yi(i) to denote the solution 
of (|2.3p in the special case where B = A. Consider the trajectories of two particles 
advected with the fluid: a reference particle with trajectory r(t), and a nearby particle 
with trajectory r{t) + dr{t). To leading order in the separation \6r\, the separation vector 
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is determined by the matrix Myi(i): we have 6r{t) — Myi(i) (5r(0). In the language of 
dynamical systems theory, a matrix with this property is termed a monodromy matrix. 
We consider volume preserving flows, so that tr[A(t)] = 0. From (|2.2p . the matrix B(t) 
also has the property that tr[B(i)] = 0, as if it were the velocity-gradient of some ficticious 
volume-preserving flow, and consequently det[M(i)] = 1. We will therefore refer to the 
solution M(i) of (|2.5p as the pseudomonodromy matrix of the flow. For the case of rod- 
like particles, where ai = 1 and a2 = in (|2.2p . it is the same as the true monodromy 
matrix of the flow. 

The degree to which the solution can be presented in closed form depends upon the 
specifics of the flow field. First we comment on the exactly solvable case of a time- 
independent flow with constant velocity gradient A, because this case already exhibits so- 
lutions showing both alignment and tumbling. In this case the matrix B is also a constant, 
and the solution of the linear auxiliary equation (|2.5p is d{t) = M(t)<io — exp(Bt)do. 
The matrix exp(Bi) may be expressed in terms of the eigenvalues and eigenvectors of B. 
The matrix B is typically non-Hermitian, and correspondingly its eigenvectors need not 
be orthogonal. The behaviour of the solution is determined by the eigenvalues A^. We 
consider incompressible flow, implying that tr[B] = 0, so the eigenvalues sum to zero. 
We describe both the three-dimensional and two-d i mensi onal cases below (the following 
discussion overlaps some comments made bv ISzeri (1993)). 



Apart from degenerate cases, the spectrum may take one of three forms in three di- 
mensions: 

(a) Eigenvalues real and distinct, with at least one of them positive. The axis of the 
particle aligns with the eigenvector m+ corresponding to the largest eigenvalue, A_|_ of B 
(in the following we refer to these as the dominant eigenvector and eigenvalue). 

(&) There may be a real and positive eigenvalue A+, and a complex pair with negative 
real part. In this case the axis also aligns with the dominant eigenvector, 

(c) There may be two complex conjugate eigenvalues with positive real part (so that 
the real eigenvalue is negative), with complex conjugate eigenvectors. When these two 
eigenvectors are combined with complex-conjugate coefficients, the resulting real vector 
lies in a plane. In the long-time limit the vector d(t) spirals outwards in this plane. This 
case corresponds to a tumbling motion of the particle. 

Another way to understand the dynamics of the vector d{t) is to write M(<;) as a 
normal form: 

M(i) = XN(i)X-i (2.7) 

where X and N(t) are real-valued matrices. In case 1, the matrix N(<) is diagonal, with 
diagonal entries exp(Aii). In cases 2 and 3, the matrix N(i) is in block-diagonal form 
with a 2 X 2 block describing a spiralling motion, 

(exp(— iAt)cos(cjt) exp(— iAt) sin(a;t) \ 
-exp(-iAt)sin(wt) exp(-|At) cos(cji) . (2.8) 

exp(Ai) / 

Here A is the real eigenvalue of B, and the complex eigenvalues are — ^A ± ioj. The 
spiralling motion may be attractive (spiralling- in, when A > 0), which is case 2, or 
repelling (spiralling-out, when A < 0), which is case 3. 

In the case of two-dimensional incompressible flow, the matrix M may have two re- 
ciprocal real eigenvalues (the hyperbolic case), or else two complex conjugate eigenvalues 
which lie on the unit circle (the elliptic case). In the hyperbolic case the vector d{t) comes 
into alignment with the eigenvector of B which corresponds to the positive eigenvalue 
(the dominant eigenvector). In the elliptic case, where B has purely imaginary eigen- 
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values zbiw, the pseudomonodromy matrix M(t) can be expressed in terms of a normal 
form, analogous to (|2.8p . with N(t) replaced by a rotation matrix R(a;t) representing 
rotation in the plane by an angle cut. 

The matrix B is typically non-Hermitian, and correspondingly its eigenvectors need 
not be orthogonal. We briefly consider the consequences of this observation in the two- 
dimensional case. If the matrix B is elliptic, the vector d{t) rotates, the linear transfor- 
mation X in (|2.7p transforms the circular motion of R(cjt)X~^<io to motion on an ellipse. 
In some circumstances this ellipse may have a large aspect ratio. In this case the vector 
n{t) = d{t)/\d{t)\ will spend most of its time nearly aligned with the long axis of the 
ellipse, reversing direction rapidly at times separated by tt/w. 



3. Order parameter and light scattering 

3.1. General definition of the order parameter 

The alignment of the particles may be described by an order parameter. If the rheoscopic 
fluid is left to stand for a while, the crystals become randomly oriented due to Brownian 
motion. When the fluid is set in motion, the crystals start to align and at later times we 
can describe the distribution of angles by a probability density. In the case we consider 
below the particles are aligned in a plane so their direction is defined by a single angle 
6. Because the direction vector is non-oriented, the probability density P{9) satisfies 
P{d + Tr) = P{9). This probability density will depend upon both position and time, but 
we suppress the arguments r and t in the discussion below. 

A suitable order parameter for the rod-like particles can be obtained from P{9) by 
first calculating the inertia tensor of the rods, which has components: 

hj^ f d9 P{9){n, ■ n{9)){n, ■ n{e)) (3.1) 
Jo 

where n{9) is a unit vector in the direction 6. The three distinct components of In, I12, 
I22 are not independent, because the vector n{9) is constrained to have unit length. They 
can be mapped to the order parameter vector C, as follows. The inertia tensor has real, 
positive eigenvalues Xi, T2 and corresponding orthonormal eigenvectors Ui, U2, with 
Ii ^ l2- The eigenvalues satisfy Xi -I-I2 — 1, and the case Ti = 1 corresponds to perfect 
alignment, whereas Ii = T2 ~ ^ corresponds to an isotropic distribution. We define 
to be a non-oriented vector in the direction Ui with magnitude which is a function of 
Xi —T2- Let us consider a special case where the rods align with the direction 9 with 
probability p, or else are randomly distributed with probability 1 — p, that is 

P{9) ^l[5{9-e)+5{9-9-7:)]+^—^. (3.2) 

Z ZTT 

It is natural to define the order parameter so that C, = pn(^) in this case. For this 
distribution, in the case ^ = we find Xi = (1 + p)/2 and X2 = (1 — p)/^, so that 
Xi — X2 = p. We therefore define the order parameter as 

C= (Xi-X2)i7i . (3.3) 

This is a general definition for the order parameter of rod-like particles in two dimensions. 
An analogous definition can be used in three dimensions, where a general inertia tensor 
has six independent components, but the inertia tensor for the rod directions has five 
parameters because of the constraint that |n| = 1. 
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Figure 5. Illustrating the geometrical construction used to determine the probability density 

for the angle, P(ff)- 

3.2. Order parameter in terms of the monodromy matrix 

Let us consider the evaluation of this order parameter for the case where the particles are 
initially randomly oriented, so that the initial direction no in (|2.2p is uniformly distributed 
about the unit circle. According to the solution presented in section [21 a vector ng on 
this circle is mapped to a vector d{t) which lies on an ellipse. This ellipse is described 
by its aspect ratio, ^ 1, and by the direction of its longest axis, 9. In the following 
we obtain the probability density P{0) and use this to obtain the order parameter in 
terms of v and 0. 

An angle interval dcj) on the unit circle is mapped to a segment of the ellipse which 
is at an angle 9 to its longer axis, and which spans an angle interval d^. The angle 6 is 
independent of the overall scale of the ellipse, and we find it convenient to consider the 
case where the short axis intersects the unit circle (see figure [5h, where ^ = so that 
the long axis is horizontal). The probability element for the direction of no being in the 
original interval is dP = d(/)/27r. This is the same as the probability element for d{t) 
being in the interval d9 on the ellipse, so that the probability density P{6) satisfies 

dP = ^d(f> = P{9)de . (3.4) 
27r 

An elementary geometrical construction can be used to surmise the relation between 
def) and d^. Instead of considering the mapping of a circle to an ellipse, let us consider the 
image of a narrow annulus of angular width dcj) between a circle with unit radius and one 
with radius 1 — e (with e <C 1), so that the area of this element is dA ~ edcj). The element 
of the annulus is the set difi^erence between two segments of discs spanned by an angle dcj), 
one of unit radius, the other of radius 1 — e. These segments are transformed into regions 
which may also be approximated by segments of circles: the larger one is approximated 
by a segment of a circle radius r spanned by an angle d^, having area ^ ^r^d9 and the 
smaller one by a segment which is smaller in area by a factor (1 — e)^ 1 — 2e (see figure 
[5)3). The area of the transformed image of the annulus is therefore dA' — er^d9. Because 
the transformation from a circular region to an ellipse stretches the x-axis by the factor 
ly, we also have dA' ~ vedcj). We conclude that d0 = r'^d9/v, where r is the distance from 
the origin to a point on the ellipse at angle 9 from the long axis. The equation of the 
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ellipse is = + iy^y^, where x = r cos 0,y = r sin6', so that = r^[(z/^ — 1) sin^ 61 + 1]. 
Using (j3.4|l we therefore conclude that the probability density for the direction of the 
vector d in (|2.4p is 

r'^ V 1 

p(e) ^ — = 5 = . (3.5) 

^ ' 2-Kv 27r (i/2 _ 1) sin2(0 - 6*) + 1 ^ ' 



Using the identities 



27r 



VA+T- 1 







da; , = 27r 

A sin"^ X + 1 A 



^'^ sin^ a; „ \/1h~T - 1 

9 = ^''^ , 

A x + 1 A\/ A + 1 



dx ^ " ^ = 2^ ^;^ ;/ / (3.6) 



we find that for this probability density the elements of the inertia tensor are 



111 = 1 - I22 = ^ cos^ Q H 5— sin^ ( 



I/+ 1 J/+ 1 

/12 = ^^-^ cos sin . (3.7) 

+ 1 

The eigenvalues of the inertia tensor are then Xi = — ^ and Xt = ^-r- The order 
parameter for an initially uniform angular distribution is therefore 

C = ^n(0-) (3.8) 

where n(0) is a unit vector in the direction d. It remains to express the aspect ratio v ^ 1 
of the ellipse in terms of the matrix M. The equation defining the unit circle |no| = 1 
can be written a; • a; = 1. In terms of x' = Ma;, this condition becomes the equation for 
an ellipse: x' ■ Ka;' = 1, with 

K = (M~i)'^M-i = (MM^)-i . (3.9) 

The aspect ratio v is therefore the square root of the ratio of the eigenvalues of the real, 
symmetric positive definite matrix K. This may also be determined from the ratio of the 
eigenvalues of — MM"''. If the matrix has eigenvalues Ai, A2 with corresponding 
orthonormal eigenvectors Ui, U2 ordered so that Ai > A2, then the parameters in 



are then v — \JXi/X2 and n(^) ~ Ui. 

In a generic flow, the matrix B(t) is neither constant nor periodic, and we expect 
that the solution of (j2.3|) will have a positive largest Lyapunov exponent. In this case, 
the pseudomonodromy matrix M will become hyperbolic almost everywhere, having a 
unique largest eigenvalue, which increases as time increases. The rods will then align 
very close to the direction of the dominant eigenvector, irrespective of their initial orien- 
tation. However, if the pseudomonodromy matrix remains elliptic, there is no dominant 
eigenvector and the final direction remains dependent upon the initial orientation. In the 
hyperbolic case where the rods approach perfect alignment, the order parameter vector 
approaches a unit vector, but in the elliptic case it is shorter than unit length. 

3.3. Relating the order parameter to light scattering 

The order parameter can be investigated experimentally by examining the reflection of 
light by the rheoscopic fluid. Because we are primarily interested in two-dimensional 
flows, we consider how the light scattering may be related to the order parameter in the 
case where the illumination is confined to a surface. By way of examples, this is relevant 
when the rheoscopic fluid is a thin layer floating on a denser, immiscible fluid, or when 
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the rheoscopic agent is used without dilution, so that the optical depth is very small 
(implying that scattered light comes from a thin layer close to the surface). The image 
contrast is greatest when the illumination comes from a direction in the same plane as 
the surface, and we choose to specify its direction by means of the angle of the direction 
perpendicular to that from which the beam is incident. 

The intensity of light reflected by the microscopic crystals depends upon their orien- 
tation relative to the direction of the source of the light. The angular dependence of the 
scattering depends upon a variety of factors, of which the ratio of the size of the crystals 
to the wavelength of light and their surface roughness are important. If the crystals are 
aligned with their long axis at angle 6, the intensity of the scattered light will be f{(j) — 9), 
for some function / which is even and periodic with period tt. In the following we consider 
the limit where the crystals are smaller than the wavelength of light, in which case the 
amplitude of the scattered radiation is proportional to the projected area of the crystal 
in the direction of the incident light. This implies that a rod at angle 9 scatters light 
from a source which is perpendicular to the direction (j) with an intensity proportional to 
cos^{9 — 0) + 7, where 7 is a contribution arising from diffuse background scattering. In 
our subsequent discussion we shall use this form for the scattering kernel, with 7 = 0. 

More detailed information about the orientation of the particles may be revealed by 
using three different light sources with different colours, illuminating the fluid from three 
different directions. The intensity of the scattering of light from a given source depends 
upon the direction of the particle relative to the direction of the light source. At any 
given position the fluid reflects with a colour C determined by additive mixing of the 
scattered light from red, green and blue (i?, G, B) sources, which we assume are arranged 
about the sample at directions separated by 120°, as illustrated in figure [H This results 
in the light being scattered with a colour C which is determined by additive mixing of 
the primary colours i?, G, B: 

C = /(0)i? + /(27r/3)G + /(47r/3)S (3.10) 



where in the limiting case of short rods I{9) is the inertia of the axial distribution relative 
to the direction 9: 



I{9)= A9' P{9')coB^{9 -9') . (3.11) 
^0 



In principle just two of the functions /(O), /(27r/2) and /(47r/3) are sufficient to determine 
the two parameters of the order parameter. However, using three colours has two advan- 
tages: with three colours the ratios of the scattered intensitie s can be used, so that the nor- 
malisation of the intensities is not relevant. Also, as shown bv lBezuglvv. Wilkinson fcMehlid 

with three colours the Poincare index of singularities can be visualised directly. 
The mapping between the order parameter C, and the colour G of the refiected light is 
illustrated in figure [D The use of coloured light sources to enhance rheoscopic images 
was previously suggested bv iThoroddsen &: Bauer (I1999I) . Their work does not consider 



the relation between the colour mixing and the ordering of the particles. 

For larger crystals the function cos^(0 — 9') is replaced by another function f{9 — 6'). 
This would make a quantitative but not a qualitative difference to the colour images 
which are displayed here. 
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4. Steady flows in two dimensions 

4.1. Hyperbolic and elliptic bands 

Now we turn to considering rheoscopic particles in a steady incompressible flow where the 
velocity vecto r is confined to a plane with coordinates (x, y). This system was previously 



considered bv lSzeril(|l993l) . who showed that there are regions (which we term hyperbolic 



bands) where the particles align with each other. His paper gives a treatment of equation 
(|2.ip using concepts from dynamical systems theory. In particular, the regions in which 
the particles align are determined by looking for stable fixed points of the Poincare map 
for the particle axis direction as it is advected around a contour of the stream function. 
Below we use the solution (|2.6|) . which gives a more thorough insight into this system, 
as well as being more computationally efficient (because it is not necessary to repeat the 
calculation for different initial directions of the rod). Together with the quadratic form 
for determining the order parameter, (|3.9p . this also allows us to describe the alignment 
of the particles in the regions outside the hyperbolic bands. 

The velocity field of a steady, incompressible two-dimensional flow may be derived from 
a stream function 7/;(a;, y): we have v = (dip/dyj—dip/dx). By analogy with Hamiltonian's 
equations of motion for a one-freedom autonomous system, we see that the trajectories 
follow contours of the stream function, so that a trajectory labelled by ipo is defined 
by writing ip{x,y) = ipo- The contours may be either closed or open. Particles which 
are advected along a closed contour have a periodic motion, with a period T (which is 
a function of tpo). This periodicity simplifies the analysis of the behaviour of advected 
particles, and we concentrate on the periodic case. (Periodic behaviour can also occur 
if ipix, y) is periodic in one or both variables, and our discussion is readily extended to 
such cases). 

We have seen that the behaviour of axisymmctric particles is determined by the pseu- 
domonodromy matrix M(i). In the two-dimensional incompressible case this matrix is a 
2x2 matrix which satisfies det[M(t)] = 1. Such a matrix is either hyperbolic, having two 
reciprocal real eigenvalues, or elliptic, with two mutually conjugate complex eigenvalues 
with modulus equal to one. The character of this matrix is readily determined from its 
trace: if |tr[M]| > 2, the matrix is hyperbolic, whereas if |tr[M]| < 2, the matrix is elliptic 
(and if |tr[M]| = 2, the matrix is a shear). By comparison with the case of constant ma- 
trix B which was discussed in section[2l we anticipate that if the matrix M is hyperbolic, 
the advected particles tend to approach a given direction, whereas the elliptic case is 
associated with tumbling motion. This expectation turns out to be correct, in a qualified 
sense as discussed below. 

We can label points on a closed trajectory by the time to taken to reach the point from 
an arbitrary reference point on the orbit. Let M(i,to) be the pseudomonodromy matrix 
for the trajectory which starts at time and at the point labelled by t^, ending a time 
t. Let us consider the evaluation of M(t,to), in the case where t is written in the form 
t = ti + NT (where T is the period and N an integer). We can express this general matrix 
in terms of a pseudomonodromy matrix for a single cycle, Mq = M(T, 0), together with 
matrices representing short time evolution for a fraction of a cycle. We can write 

M{t,to) = M(ti,0)[Mo]^M-i(to,0) . (4.1) 

This shows that the long-time behaviour is determined by the character of the matrix 
Mo, which can be computed by propagating a solution of (|2.5p for a finite time. In 
particular, if Mo is hyperbolic, the matrix M(t,to) will have one eigenvalue which is 
much larger than the other when t — to — > oo. Because the eigenvalues of a matrix are 
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invariant under a similarity transform, the structure of (|4.ip implies that the character of 
Mo (hyperbolic or elliptic) is independent of the choice of starting point on the contour. 

Because the elliptic or hyperbolic character of a trajectory is independent of its starting 
point, we can label the contours of the stream function according to the character of the 
one-period monodromy matrix, Mq. From (|4.ip . we see that when ti = and when Mq 
is hyperbolic, the particles align with the eigenvector of Mq corresponding to its largest 
eigenvector. More generally, in the hyperbolic case the orientation at any position aligns 
with the direction of the dominant eigenvector of the monodromy matrix for the one- 
period orbit which ends at that position. On contours where Mq is elliptic, at any given 
position the particles continue to tumble as t — to —^ oo. The contours of ip{x, y) are may 
therefore be divided into elliptic bands, where the pseudo-monodromy matrix is elliptic 
and the particles tumble, and hyperbolic bands, where the pseudo-monodromy matrix 
is hyperbolic and where the direction approaches a constant vector field. These bands 
are analogous to the bands wh ich occur for the solution of the Schrodinger equation 
for a one-dimensional potential ( Ziman 19761 ). where the transfer matrix for solutions of 
the Schrodinger equation plays the same role as the monodromy matrix for a single orbit 
Mq. The hyperbolic bands correspond to the band-gaps in the solution of the Schrodinger 
equation, where the wavefunction increases exponentially in one direction, so that there 
are no satisfactory eigenstates. The elliptic bands correspond to the energy bands of the 
Schrodinger equation, where its solutions are Bloch waves ((ziman 1976). We remark that 
our discussion may also be viewed as an example of the application of Floquet t heory. 
We em phasise that the hyperbolic bands are the same as the aligning regions in ISzeril 
(1993). 

We remark that in the special case where the particles are rod-like (that is, q;2 ^ in 
(j2.2p ). the matrix M(i,to) is the true monodromy matrix. The one-period monodromy 
matrix for a periodic two-dimensional flow is always a simple shear, and rod-like particles 
will always align with the contours of the stream function. 

The results in section [2] above show that in a simple shear flow, the particles always 
tumble rather than coming into alignment (except for the limiting case where the aspect 
ratio of the rods is infinite). Because a steady two-dimensional flow locally resembles a 
shear flow, it might therefore be expected that the transfer matrix Mq would always be 
elliptic, because it can be thought of as a product of matrices each of which would indi- 
vidually be generated by an elliptic flow. This need not be the case, however. It is known 
from studies of Anderson localisation for the one-dimensional Schrodinger equation that 
products of elliptic matrices can be hyperbolic (Mott & Twose 1961). (In the context 
of Anderson localisation, this statement is equivalent to the observa tion that loc alised 
states occur for energies where there is no classical potential barrier ( Zimaniri976f )). We 
therefore conclude that alignment of particles around periodic trajectories is possible, 
although it might be argued to be un-expected. 

We now turn to consider an example of textures formed by the alignment of axisym- 
metric particles in steady two-dimensional flows. Figure [5^ displays the contours of the 
stream function for a two-dimensional flow, exhibiting saddle points as well as extrema. 
The flow is a 'journal bearing' flow, where the circular boundaries rotate with different 
angula r veloc it ies. The s tream function for this flow was obtained by Jefferv ( 19 22al ). 

fTln 



Miilleil (Il942[ ). IWanilie^ (|l95nl ). and is discussed in detail in (jBaUal fc Rivlinlll976l )7ln 



this example, the walls rotate in the clockwise sense, with the angular velocity of the 
inner wall exceeding that of the outer wall by a factor of 20. The radius of the in ner wall 
is 0.3 times that of the outer wall, and the eccentricity parameter e of (jBallal fc Rivliiil 
1976h is |, so that the centre of the inner boundary is offset by a multiple of 0.525.. 



times the radius of the outer wall. This system can be realised physically by filling the 




Figure 6. a Contours of the stream function tp{x,y) for a journal bearing system: both circular 
boundaries rotate in the same direction, with the angular speed of the inner boundary exceeding 
that of the outer boundary by a factor of 20. b These contours can be coloured according to 
whether the transfer matrix Mo is elliptic |tr(Mo)| < 2 (red), or hyperbolic, tr(Mo) > 2 
(blue) and tr(Mo) < —2 (green). Each hyperbolic band is labelled by its Poincare index. In 
this illustration we set qi = 0.95, Q2 = 0.05 in equations (12. ip . (|2.2|l . (which corresponds to 
ellipsoidal particles with aspect ratio f3 = vTQ — 4.36..). 



space between two vertical rotating cylinders with a rheoscopic fluid, and figure |4] (which 
will be explained fully in section is an illustration of the complexity of the pattern 
of light scattering from the surface which could be observed in the long-time limit. The 
hyperbolic bands (blue or green) and elliptic bands (red) are illustrated in figure[6jD, with 
the hyperbolic bands shaded blue if tr(Mo) > 2, green if tr(Mo) < —2 (the reason for 
making the distinction between thee two hyperbolic cases will be considered in section 
15. 3p . There is a contour which marks a transition from trMo > 2 to trMo < 2 with- 
out passing through an elliptic zone: this is possible because the contour is a separartix 
where the topology of the contours chaii g es. Th is example of a steady two-dimensional 
flow is the same as was studied by ISzeril (|l993h . and the hyperbolic bands in figure [6Jd 



correspond to the aligning regions which were obtained by Szeri. The Poincare indices of 
the hyperbolic bands ar e also s hown , and these are equal to one-half of the 'flip numbers' 



which were discussed in ISzeril (|1993( ) 



4.2. The order parameter in elliptic bands 

Let us consider the form of the order parameter in the elliptic bands. Provided the period 
T of an orbit depends upon the stream function a passive scalar function will be wound 
into an increasingly tight spiral under the action of a two-dimensional steady flow. Its 
lines of constant density will become closely aligned with the contours of the stream 
function, with the scalar having an approximately periodic behaviour when traced in a 
direction perpendicular to the streamlines. Figure [1] showed an example of the evolution 
of the order parameter as time increases, showing the development of an increasingly 
tight spiral pattern. However, we shall see that the behaviour of the order parameter 
is more complicated than that of a passive scalar, in that its variation in a direction 
perpendicular to the contours of ipix, y) is quasiperiodic rather than periodic. 

Consider the variation of the order parameter within an elliptic band as a function of 
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\^ NT(i,o) = t 

(N + l)T(i,o- Ma)=t 

Figure 7. Illustrating the coordinates Ai/), r which are used in the discussion of elliptic bands. 



position for large time t, in the vicinity of a reference point which lies on a closed contour 
of V'- In the neighbourhood of this reference point (a;o,?/o) we use two coordinates I^ip 
and T to label points {x,y). We define Aip = ip{x,y) — tp{xo,yo)- We define a reference 
point on other contours of "0 by drawing a hne which is perpendicular to the contour 
passing through (xo,yo)- We label the distance along a contour by the time r taken to 
reach that point starting from the reference point on the orbit. This coordinate system 
is illustrated in figure [T] 

Now let us specialise by taking the reference point to lie on a contour such that t is a 
multiple of the period T, so that t = NT for some integer N. For a set of isolated contours 
the motion will also be periodic, making a different number of orbits in the same time t. 
For large t these contours are approximately evenly spaced, with the spacing A-^o of the 
contours of the stream function being 







1 A' 


TVdT 




NA" 



(4.2) 



where A{ip) is the area enclosed by the contour with stream function ip. 

The transfer matrix may be written in terms of its normal form, similar to (|2.7[) . First 
consider the form of this matrix along the line r = 0. At (x, = {xo,yo), we have 
M(0, 0) = M^, where Mq is the transfer matrix (that is, the pseudomonodromy matrix 
for one orbit). We write the transfer matrix in normal form as follows: 



Mo = XR(6'o)X- 



(4.3) 



where R(6') is a rotation matrix for angle 6. When ip changes by A-^o, the trajectory 
makes one additional orbit, so that the transfer matrix becomes M(A?/'o,0) — M^'*'^. 
We can therefore write M{Aip,0) = X R(6') X^^ Z(A?/'/A?Ao), where Z{x) is a 2 x 2 
matrix which is a periodic function of x, with 



= (70 



A^ 



and 



Z{x + l) = Z{x) , Z(0) = I. 
With these notations and definitions, for a general position the transfer matrix is 

M(Ai/;,t) = M(r)XR(6')X-i Z(AV'/AV'o) M-\t) . 



(4.4) 



(4.5) 



(4.6) 



In the limit as — > oo the order parameter depends increasingly sensitively upon ^, 
but the sensitivity to r is independent of A^. The dependence of upon Aip is quasiperi- 
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Figure 8. Illustrating the evolution of the angular variables in the representation (|4.6|l . 9 and 
At/j, as the contour label tpo is varied. The periods of these variables are 2n and At/jq respectively, 
and (|4.4p implies that the slope of the line is 9o/Aipo. The evolution can be 'folded' into a unit 
cell, and provided Oq /2n is an irrational number this reduced dynamics is ergodic. 

odic, being associated with two periods. One period Aipo is associated with the change 
in required for the trajectory to make an additional orbit in time t. There is another 
periodicity associated with the change in required for the phase 6 in (|4.4p to increment 
by 27r. This additional periodicity is AV'i — 27rA-0o/^o- 

Equation (|4.6p can be used to obtain a complicated expression for the matrix of the 
quadratic form describing the order parameter, l|3.9p . which is arguably too unwieldy 
to be of much use. However, in the next section we shall see that although the order 
parameter depends increasingly sensitively on position in the long-time limit, the order 
parameter of the locally-averaged orientation has a very simple representation. 

5. Averaging, singularities and topology of the order parameter 



We have seen that in the elliptic bands the order parameter varies increasingly rapidly 
as a function of ^p in the limit as t oo. Eventually the order parameter fluctuates 
on a length scale which is small compared to the resolving power of the eye. In this 
limit it is necessary to perform a local average of the inertia tensor (j3.ip representing 
the distribution of orientations. The order parameter of this locally averaged quantity 
determines the appearance of the rheoscopic suspension in the long-time limit. 

In the limit as t ^ oo the periods associated with varying ip^ namely A-^i and AipQ 
respectively, both approach zero. As the contour -00 is varied, the values of 9 and Aip 
both change linearly, along a trajectory illustrated in figure |51 The local averaging of 
the orientation distribution is effected by averaging along this trajectory. Because of 
the periodicity, the trajectory can be 'folded back' into a single unit cell. The folded 
trajectory will fill this unit cell provided 9^/21^ is an irrational number (that is, not a 
ratio of two integers). Because rational numbers are a measure zero case, we may perform 
the local average by averaging (|4.6p over the unit cell in figure [51 

Consider the behaviour of the order parameter of the locally averaged orientation in 
terms of the representation (|4.6p (without loss of generality we may consider the line 
T = 0). We consider a region which is large compared to both of the periods A-0o and 



5.1. Local average of the order parameter 
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n X-^Zn (R(e)X-^Zn)9 X(R(e)X-^Zn)fl 

Figure 9. Illustrating the transformations which are applied in succession to the distribution 
of the initial direction vector n, in order to produce the vector d = M n, where M is expressed 
in the form (|4.6p . The vector n is initially randomly distributed around a unit circle (a). After 
application of the transformation Z, this circle is transformed into an ellipse, with the 
parameters of the ellipse depending periodically upon Atp (b). If we average over the rotation 
angle of the matrix R(S), the vectors which are randomly distributed on an elliptical curve are 
mapped into an annular region (c). This region is mapped into an elliptic annulus by the final 
transformation X (d). 



Aipi (note that both periods approach zero in the long-time limit, so this region can be 
made arbitrarily small). The orientation of n is initially distributed randomly around the 
unit circle. The matrix Z{Aip / A-ipo) maps this circle to an ellipse, the parameters of 
which depend periodically upon Aip, with period Atpo (this is illustrated schematically 
in figure [S^,b). We will average over the period Aipo as the final stage of our argument. 
This ellipse is rotated by the angle 0, which depends increasingly sensitively on ijj in 
the long-time limit, with a period A-^i which is inversely proportional to time, so that 
we can average over the rotation angle 9. Upon averaging over 0, the ellipse is therefore 
transformed into a circularly symmetric distribution in the plane, as illustrated in figure 
The action of the matrix X transforms this annular region into a region bounded 
by two similar ellipses; see figure [5)1. These have an aspect ratio v which is the square 
root of the ratio of the eigenvalues of XX""", as described in section [31 The arguments 
developed in section [3] show that the angular distribution P{9) depends only upon the 
aspect ratio of the ellipse, and not upon its overall scale. Furthermore, although the 
radial distribution in the circular region depends upon Atjj, it is only the aspect ratio of 
the elliptic region which matters, and this is determined solely by the matrix XX'^, so 
that the average over A^ is trivial. 

Thus we conclude that in the elliptic regions the locally averaged order parameter 
approaches a limit which varies smoothly as a function of position. The averaged order 
parameter (C)('') is determined by the matrix X(r) which occurs in the definition of the 
normal form (|4.3p . in the same manner as the un- averaged order parameter C(''j^) is 
determined from the pseudomonodromy matrix M(r,t, io)- In particular, the equation 
(|3.9p for the matrix K defining the quadratic form for the inertia tensor of the angle 
distribution is replaced by 

K-^=X'X.^. (5.1) 
The locally-averaged order parameter ((J) points in the direction of the eigenvector cor- 




Figure 10. a Illustrating the locally averaged order parameter field computed using (|4.3p and 
(|5.1[1 . The hyperbolic and elliptic bands are separated by red lines, and the positions of zeros 
of the order parameter field are indicated by green dots for zeros with Poincare index equal to 
|, green crosses for zeros with index — |. b Is the same image as figure |4l with the boundaries 
between hyperbolic and elliptic bands indicated by solid black lines (and the positions of zeros 
of {^) are also marked). 

responding to the largest eigenvalue of XX'^, and if the square root of the ratio of 
eigenvalues of this matrix is fi, then \ \ — {fi ~ 1) / {fi + 1) . Equation (|5.ip is one of the 
principal results of this paper, since it expresses the long-time limit of the alignments of 
particles in terms of the normal- form of the transfer matrix Mq. The locally averaged 
order parameter field is illustrated in figure [TOb for the same journal bearing example as 
figures [4] and [6l 

5.2. Continuity of the averaged order parameter 

In a hyperbolic band, where the particles approach a fixed alignment, the asymptotic rod 
direction at any point r approaches the dominant eigenvector u+ of the transfer matrix 
Mo(r), for a periodic orbit which ends at r. In the long-time limit, a local average of 
this order parameter field, (C)('")j is also a smooth function of position throughout the 
elliptic region. We should consider whether the locally averaged order parameter varies 
conti nuously upon passing between e l liptic and hyperbolic regions. 



In IWilkinson. Bezuglvv fc Mehligj (|2008r ). we showed that the transfer matrix at a 
boundary between elliptic and hyperbolic regions where tr(M) = 2 is in the form of a 
generalised shear: 

M = R(0)S(k)R(-</.) (5.2) 
where R((^) is a rotation matrix and S(k) is a shear of the form 

S(^) = J iy (5.3) 

It follows that eigenvectors of the monodromy matrix become co-linear as we approach 
the boundary between elliptic and hyperbolic regions: both eigenvectors approach u = 
(cos 0, sin 0), while both eigenvalues approach unity. As we approach such a boundary 
from the hyperbolic side, the order parameter field aligns with this common eigenvector. It 
will prove useful to express (|5.2p in component form: introducing the notations c = cos (j), 
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s = siiKji, we find that 

KCS 



^ y —Ks"^ 1 — KCS j ' i^-^) 

However, it is not immediately clear what happens as we approach the boundary 
from the elliptic side, where the transfer matrix can be expressed in the form ()4.3|) . As 
the boundary is approached, the angle 6*0 in (|4.3p approaches zero, because tr(M) = 
2cos0o ^ 2, and in the following discussion we treat a-s a small number. It is clear 
that the matrix X in the representation (j4.3p must become singular in order to approach 
(|5.2p as 6*0 0. Let us assume that in this limit X takes the form: 



cos(0 + S(f>) cos{(j) — S(j)) \ _ f c — S(f) s c + 5(f)s 
sm{(j) + 5(f)) sin((/) — Scj)) J \ s + 5(f)c s — Scjic 



0{5r) (5.5) 



where we use the notations c = cos((^), s — sin((^), and where we shall assume that the 
small change in the angle is 

6c^=^ + 0{el). (5.6) 

We find det(X) — 29o, so the assumed form for X does indeed become singular as 0o ^ 0. 
Inserting the ansatz ()5.5|) . (15. 6p into ()4.3p . approximating 

R(0o)=I + M + O(02) , \ J) (5.7) 



and ignoring 0{9^) terms, we find: 



In the limit as we find that this expression agrees with (|5.4p . which confirms that 
the ansatz ()5.5p . (|5.6p was correct. We can now use this expression for X in equation 
(|5.ip to calculate the form of the matrix defining the quadratic form characterising 



the order parameter: we obtain 



K-^ = 2(^;:^ ^2j+O(0o')- (5.9) 

The term which is independent of 0o is a singular matrix: its eigenvectors are Ui = 
(cos 0, sin (/)) with eigenvalue Ai = 1, and U2 = (sin 0, — cos ^) with eigenvalue A2 = 0. 
This shows that in the limit as 6*0 —> the ellipse which is defined by the quadratic 
form (XX'^)^^ degenerates into a line, which is aligned with the common eigenvector of 
(|5.2p . Because the aspect ratio of the ellipse approaches infinity, the modulus of the order 
parameter approaches unity as the boundary is approached. We conclude that the locally 
averaged order parameter is continuous at the boundary between elliptic and hyperbolic 
regions (although it clearly has discontinuous derivatives). 

Finally we comment on the nature of the discontinuity of the order parameter at the 
boundary between the elliptic and hyperbolic bands. The fact that the 0(6'o) term in (|5.9p 
is equal to zero implies that the determinant of XX""" is det(K~^) = O(0§), implying 
that A2 = O(0q). This implies that the aspect ratio of the ellipse is ~ 9^'^ . Because 
tr(M) = 2 cos 6*0 has a linear dependence upon the distance d from the boundary with the 
hyperbolic region, we conclude that ''^ Vd, so that v ^ 1/d. This in turn implies that 
the magnitude of the order parameter approaches unity linearly upon approaching the 
boundary of an elliptic band, implying that {C){r) has a discontinuous first derivative. 
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5.3. Poincare indices of the order parameter 

The Poincare index must be the same for any curve lying in a hyperbohc band , be- 
cause u+ cannot have any singularities there ( Wilkinson. Bezuglyy fc Mehlid 20081 ). The 
Poincare index is most efficiently determined by evaluating it+ around a given contour 
of V' within the hyperbolic band. These Poincare indices are indicated for each of the hy- 
perbolic bands in figureElb. They were also evaluated bv lSzeri (1993): his 'flip numbers' 
are twice the Poincare index. 

We have seen that each hyperbolic band is associated with a Poincare index, and 
simulations confirm that the Poincare indices of different hyperbolic bands need not be 
equal. We have also seen that (^) is continuous everywhere, so that a Poincare index can 
also be ascribed to the averaged order parameter field (C)('") in the elliptic bands. This 
raises the following question: is there a rule for determining the difference between the 
Poincare indices of the hyperbolic bands in terms of a property of the intervening elliptic 
band? 

The analogy with Bloch bands in solid-state physics suggests that a rule for Poincare 
indices might be found. The wavefunction of a Bloch band at any given energy is char- 
acterised by a Bloch wavevector k, such that on traversing one period L of the potential 
the wavefunction accumulates a phase factor exp(iA:L). The phase 9o in ()4.3p corresponds 
to kL in the Bloch wavefunction. The wavevector is related to the monodromy matrix 
by |trM| = 2cos(A:L). On traversing a band, the wavefunction therefore rotates by tt for 
every period of the potential. By analogy, in a steady flow we might expect that the axis 
rotates by ±7r on crossing every elliptic band, which would imply that the Poincare index 
changes by ±^ on crossing every elliptic band. The following argument shows that this 
physical intuition is partially correct. 

A rule for changes of the Poincare index is obtained by the following argument. We 
assume that the elliptic region has a simple annular topology, although cases where an 
elliptic region has two or more 'holes' occur. In the elliptic band the eigenvalues of Mq 
are complex numbers, with both the eigenvalues and eigenvectors occurring as complex 
conjugate pairs. We can multiply the eigenvectors by complex numbers chosen so that 
these vectors are purely real at the band edges. Let us combine the two eigenvectors 
Ui, U2 = ul to yield a real- valued vector a = ^[iti + 1*2]. This vector depends upon 
position, because the matrix Mq depends upon the position r. In in the following we 
use the same coordinates as in section |4j so that positions within the band are labelled 
by the value of the stream function contour, -00: and by the time r taken to reach 
the point from a specified starting point on the contour (see figure [7|) . Note that the 
matrix Mq at different points around the contour is related by a similarity transfor- 
mation: Mo(V'o, t) = ^{iPott) Mo(V'o, 0) M~-^('i/'o, t), implying that eigenvectors satisfy 
Ui{'4'o,T) = M('0o, T)Mi(-0o, 0). Now let us consider some properties of the vector field 

A(0o,r) =M(^o,r)a(Vo,0) = iM(0o,r)[Mi(Vo,O)+tt2(V'o,O)] . (5.10) 

We note the following properties of this vector field: 

(a) At the inner and outer edges of the elliptic band (we label these contours 4'i and 
V'2 respectively), the two eigenvectors Mi, U2 become colinear, and the real- valued vector 
A{'ip, t) corresponds to the single eigenvector of the monodromy matrix. The vector A 
therefore corresponds to the long-time limit of the order parameter at the inner and 
outer edges of the elliptic band, and the Poincare index of A on the inner and outer 
edges corresponds to the Poincare indices (iVi and N2 respectively) of the surrounding 
hyperbolic bands. 

(6) The vector field ^(-00, r) is clearly a smooth function of position within the elliptic 
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band. Also, because M(?/'o,t) is non-singular, and the vector a(?/'o,0) does not vanish 
for any value of -00 in the interval [ipi, tp2], this vector field a{ipo, t) has no zeros in the 
elliptic band. 

(c) Let us consider a closed curve which is composed of the line r = traversed from 
the outer edge to the inner (from ipQ = 'ip2 to ipo = ipi), the inner edge of the elliptic 
band (that is, the line = V'l) traversed clockwise around one period, the line t — 
traversed from ipQ — ^pi to the outer edge "00 — find then the outer edge (the line 
tpQ = 1P2) traversed counterclockwise back to the starting point. This path is illustrated 
in figure [TT] Because the vector field A has no zeros and is everywhere smooth within 
this region, the Poincarc index N of this field evaluated on the specified path is equal to 
zero. 

(d) However, we note that the vector field A{ip,T) is periodic on the segments which 
correspond to the inner and outer edges of the elliptic band {ipo — tpi ot ipo = ^2), so that 
we can talk about a Poincare index defined on these segments of the path in isolation. 
Furthermore, because A corresponds to the order parameter field (^) at the band edges, 
we see that the contribution to the Poincare index of A on the closed paths which 
arise from the inner and outer band edges is equal to the difference between the Poincare 
index of the order parameter field at the inner and outer edges of the elliptic band. The 
can therefore deduce this difference (that is, iVi — N2) by evaluating the contribution to 
the Poincare index which arises from the two segments along the line r = 0. 

(e) Although the path in space which is followed by the two 'radial' segments of path 
considered in 3 above is the same, the vector field differs because in one case the matrix 
M(-(/;o, T) — Mo (■00, 0) has been applied to the vector a(0o, 0). This vector is constructed 
from the two complex-conjugate eigenvectors of Mq, for which the corresponding eigen- 
values may be written as exp(iX), where tr(Mo) = 2 cos{K). If we write the eigenvectors 
of Mq in the form u = a + ib, where b is a real-valued vector, then we can express the 
relation between the vector A on the two radial components of the closed path as follows: 

^(V'o, T(0o)) = cos(if (0o))a(V'o, 0) + sin(if (0o))b(V'o, 0) . (5.11) 

(f) Equation (|5.1ip leads to two possible conclusions. The band edges correspond to 
points at which tr(Mo) — ±2 — 2cos(i<'). This implies that sm{K) = at the band 
edges and cos(Ar) = ±1. If tr(Mo) has opposite signs at the two band edges, then 
equation (|5.1ip implies that A changes sign when the two radial elements of the closed 
path in figure [TT] are traversed in opposite directions. Because the Poincare index for the 
composite path is equal to zero, this change of sign implies that the Poincare indices of 
the inner and outer band edges differ by ±i. Conversely, if the sign of tr(Mo) is the same 
on the inner and outer band edges, then the Poincare indices of the inner and outer bands 
edges are equal. A more detailed argument, requiring information about eigenvectors of 
Mq as well as its trace, is required to establish the sign of the change in the Poincare 
index. 

In the solid-state physics context, the structure of the Schrondinger equation implies 
that the trace of the monodromy matrix always does change sign upon crossing a band. 
In the problem we consider here, tr(Mo) need not change sign upon crossing an elliptic 
band, so that the surmise about the Poincare index based on the solid-state physics 
analogy is only partially correct. 

5.4. Singularities of the order parameter 

As pointed out in sections [5?2l and [531 above . the locally averaged order parameter vector 
(C) varies smoothly and is defined everywhere within the elliptic bands. However we have 
seen that the Poincare index of the order parameter field may differ by ± ^ between the 




inner and outer edges of the band. When these Poincare indices are different, there must 
be at least one singular point inside the band, where there is a zero of the averaged order 
parameter vector field {C){r). We now consider the structure of the se singularities. A 
similar argument is presented in Bezuglvv. Wilkinson fcMehlig ( 20091 ). where we discuss 



singularities of the order parameter ^ for random flows. Here we discuss singularities of 
(^) for steady flows, and find that the mathematical structure of the singularities is the 
same, although the argument has a different structure. 

In the case where X in equation ()4.6p is a unit matrix, there is a singularity where 
the orientation remains uniformly distributed, implying that the order parameter vector 
vanishes. We now examine the structure of the position-dependence of this matrix in the 
vicinity of this singular point. A general 2x2 matrix A can be written in the form 

A = aR(0)diag(A,A"i)S(K) (5.12) 

described by four parameters a, 0, A, k, where S(k) is the shear matrix, (|5.3p . Consider 
the use of the representation (|5.12p to parametrise the matrix X in (|4.3p . First note 
that because the scaling constant a and the rotation matrix R(0) both commute with 
R(6'), if we express X in the form (|5.12p . the values of a and (j) are irrelevant, so that 
we may write X as a member of a two-parameter family: X = diag(A, A~^)S(t). By a 
linear transformation T of the coordinate system, we may represent the position r in the 
vicinity of a zero at ro in terms of coordinates X = {X, Y), writing X = T(i — ro). This 
change of coordinates is non-inverting (that is, det(T) > 0) and is determined so that 
A = 1 + + 0{X'^), K = sY + 0{X'^), with the sign s = ±1 chosen so that det(T) > 0. 
The position dependence of the matrix X may therefore be parametrised as 



F sY 
1-iX 



0{X^) . (5.13) 



The parameter dependence of the matrix K ^ = X X""" is therefore of the form 

K-i=(\+/ ^'^^)+0{X'). (5.14) 



This matrix has eigenvalues A± — 1 ± i?, where R = y^X"^ + Y'^, and if we write (X, Y) = 
(i?cos0, i?sin8), we find that the eigenvector corresponding to the largest eigenvalue, 
1 + R, has angle 9 = s^Q. The aspect ratio is i/ = + R)/{1 - R) = 1 + R + 0{R'^). 
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Figure 12. Illustrating the normal forms for the zeros of the locally averaged order parameter 
field {(^){X, Y): a s — +1 leads to a core singularity, b = — 1 leads to a delta singularity. 



The magnitude of the order parameter is then \C\ = R/2 + 0{R^), so that the locally 
averaged order parameter is 

(0(X,r) = |n(|e) + O(X2). (5.15) 

The field {C){X,Y) is illustrated in figure [12] for both choices of the sign s. In both 
cases the normal form of the singularity, (15.151). resem bles forms which are seen in ridge 



patterns of fingerprints (first described bv lHenrv ( 19001 )): we have a core singularity when 
s = +1 or a delta singularity when s = — 1. 

The singularities of our order parameter field are very closely related to 'umbilic points' 
on surfaces, where the height z above the Cartesian plane is z = f{xi,X2)- An umbilic 
point is a point where the magnitudes of the principal curvatures are equal, so that the 
su rface is locally isotropic . Different ways of categorising umbilic points are discussed 
bv lBerrv fc Hannavl (|l977l ). The principal curvatures are the defined by the eigenvalues 



and eigenvectors of the real symmetric Hessian matrix, with elements d'^f/dxidxj. This 
is analogous to considering the matrix K(a;,?;) discussed above, and the direction of 
one of the principal axes of curvature has the same singularities as the direction of 
our order parameter. T he classification of directions of principal curvatures discussed by 
Berry &: Hannavl (|l977 ) lists three types of singularity, star, lemon and monstar. The star 



is equivalent to the delta singularity of fingerprint patterns. The lemon and monstar are 
subdivisions of the core singularity. They are distinguished by the number of lines along 
which the vector field is aligned radially. In the core and delta singularities illustrated 
in figure [T^l there is one such line for the core singularity and there are three such 
lines for the delta singularity. This figure illustrates the two singularities expressed in 
the normal form coordinates, {X,Y). Upon transforming back to the original Cartesian 
coordinates, x = T^^ X, however, angles need to be preserved, and for some choices of 
T the core singularity has two additional lines where the vector field po i nts ra dially. Core 
singularities with one radial line are termed le mons i n Berry fc Hannav ( 19771 ). and those 



with three radial lines are monstars. lOennisI ( 2008 ) gives a clear and nicely illustrated 
discussion of monstar singularities. 

Zeros of the order parameter can be identified in the journal bearing example and their 
positions are plotted in figure [TOl These are generic zeros with the same structure as the 
normal forms discussed above, but the transformation T from the original coordinate 
system to that of the normal forms is close to being singular, so that the structure of the 
normal forms is highly distorted in figure [TOl The particles in contact with the moving 
walls tumble, but their order parameter is in alignment with the walls at the boundary. 
The Poincare index of both the inner and outer boundaries is therefore +1, like that of 
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the 'vortex' in figure [2]a. The Poincarc index of the hyperbohc regions can be seen to 
obey the rule discussed in section r5.3( changing by no more than ±^ on crossing each 
annular elliptic band. There are a total of sixteen zeros of the order parameter lying in 
the elliptic bands: eight cores and eight deltas. Their topological charges can be seen 
to be consistent with the changes of the Poincare index on crossing elliptic bands. It is 
interesting to note that in figure [TOl the set of zeros is symmetric under reflection, despite 
the fact that the order parameter field is not. We discuss the behaviour under reflection 
in section [521 below. 

We also investigated the alignment of particles for a 'generic' stream function, defined 
by a two-dimensional real-valued Fourier series on a square domain with random Fourier 
coefficients. An example is shown in figure fT3l which shows contours of the stream func- 
tion (a) , regions where the transfer matrix is hyperbolic and elliptic and Poincare indices 
of the hyperbolic bands (b), the averaged order parameter field (c) and its colour map- 
ping (d), with the zeros marked. It can be verified that this example satisfies the rule 
discussed in section [5731 which related the Poincare indices to tr(Mo). 

5.5. Behaviour close to centres of rotation 

The arguments in section [574l above show how the existence of zeros of the order parameter 
may be deduced in elliptic bands which lie between hyperbolic bands. We now discuss 
the regions surrounding stable fixed points of the fluid flow. 

The fluid flow has elliptic fixed points (centres of rotation) at maxima and minima of 
the stream f unction ip{x, y). Szeri analysed the motion of rods at these elliptic fixed points 



(|Szerilll993[ ). He showed that the rod axis rotates at the fixed point, with a frequency 
which is less than the frequency at which fluid elements rotate at this point. This can 
be seen immediately from the general solution ()2.6p . because at the fixed point the 
pseudomonodromy matrix is constant in time, so that the solution discussed in section 
[Hcan be applied directly. It follows that the elliptic fixed points of the flow are always 
surrounded by elliptic bands of the pseudomonodromy matrix. 

It is natural to ask whether the stable flxed points of the fluid flow correspond to zeros 
of the order parameter. At the stable fixed point, the pseudomonodromy matrix M is 
generated by exponentiating a constant velocity gradient A, so that the transfer matrix 
is Mo = exp(Br), where T is limit of the period of the flow fluid as the fixed point 
is approached, and B is the matrix defined by (|2.2p . evaluated at the fixed point. The 
normal-form decomposition of Mq will be a pure rotation if the minimum or maximum 
is (to leading order) circularly symmetric, but in the general the matrix X which occurs 
in ()4.3p will not be the identity matrix. Thus we see that, except where fixed points are 
isotropic, the order parameter is non-zero at stable fixed points of the velocity field. 

Inspection of figures [TOl and [131 confirms that the stable fixed points always occur in 
regions where the transfer matrix is elliptic, and that stable fixed points do not coincide 
with zeros of the order parameter. 

5.6. Reflection symmetry 

The journal bearing example which is illustrated in figures HI [6l and [TOl above has a 
stream function which is invariant under reflection about the line y = 0. It is interesting 
to consider the extent to which this symmetry is reflected in the alignment of the rod-like 
particles. This is most easily understood by comparing the transfer matrix Mq at a point 
r — {x,y) with its value Mg^ at a reflected point = {x, ~y)- We find it convenient to 
represent the effect of the reflection by a matrix S: 

r-^^Sr, S=( J _M . (5.16) 
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Figure 13. Illustrating particle alignment in a randomly-generated stream function: a Contours 
of the stream function, which is periodic on a square of length i. b Shows the elliptic bands, 
|tr(Mo)| < 2 (red), and hyperbolic bands, tr(Mo) < -2 (green) and tr(Mo) > 2 (blue). The 
hyperbolic bands are labelled with their Poincare index, c Shows the locally-averaged order 
parameter field in the long-time limit, d Shows the colour mapping of the locally-averaged 
long-time order parameter field. The zeros are marked with green dots (cores) and crosses 
(deltas). In this figure the aspect ratio parameters of the rod- like particles are ai = 0.875, 
02 = 0.125, (3 = V7- 



The sense of rotation (clockwise or counter-clockwise) about a contour of the stream 
function is reversed under reflection, which corresponds to taking the inverse of the 
pseudomonodromy matrix. We can therefore construct Mg^ by applying a reflection, 
applying time- reversed propagation at r, and then reflecting again, that is 

M^ = SMoS. (5.17) 

In component form, the elements two transfer matrices are therefore related as follows: 



Mo 



mil rni2 

^21 m22 



TO22 mi2 

TO21 mil 



(5.18) 
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The corresponding matrices describing the quadratic form for the time-averaged order 
parameter, K^^ — XX"^ at r and (K^)~^ at are therefore related as foUows: 

K-i = ( ^2 A (j^R)-i ^ ( fcn -fci2 \ (5^^9) 

Since they have the same determinant and trace, they have the same eigenvalues and the 
length of the order parameter vector is the same at the reflected point. The direction of 
the order parameter at two points satisfies 6 + 0^ = 27r, because the signs of off-diagonal 
elements are opposite. This relation between the directions at reflected points is apparent 
in flgure [TOb.. 

It also follows that zeros of the order parameter come in symmetric pairs, except 
where there is a zero on the axis of symmetry. Also, the Poincare index of a zero and 
its reflected partner must be the same. This is conflrmed by inspection of flgure [TUb . A 
further consequence is that the directions on the symmetry axis can only be aligned with 
the axis, or else perpendicular. 



6. Concluding remarks 

The alignment of small anisotropic particles due to velocity gradients in fluid flows is 
a signiflcant problem, with a broad range of potential applications. There is as yet no 
general solution for a triaxial body in a three-dimensional flow, but the special case of an 
axisymmetric body is expected to exhibit most of the physically important phenomena. A 
simple and powerful general solution for the axisymmetric case was given by Szcri (1991, 
who showed how the orientation may be obtained from a companion linear problem. This 
present work is the third of three papers which have investigated the consequences of 
this solution in different situations. These concluding remarks will set the results of this 
paper in context with our earlier work. 

The characteristics of the solution depend upon whether the flow is chaotic or recircu- 
la ting, and upon whether we average ove r a random initial conflguration of the particles. 
In I Wilkinson. Bezuglvv fc we considered the case where the initial orien- 

tation of the particles is not random, and where they are advected in a non-steady flow 
(which may assumed to have a positive Lyapunov exponent in most cases). We showed 
that the particle orientation fleld n(r,i) is, strictly speaking, a smooth function of the 
position r, but that numerical simulations in two dimensions exhibit apparent singulari- 
ties, which resemble the core and delta singularities in the ridge patterns of flngerprints. 
We showed how the occurrence of these apparent singularities can be explained. We also 
discussed the behaviour of the solution (|2.6p in the long-time limit: we showed that, de- 
spite the increasing sensitivity of M(r, ro, t) to the flnal position r, the direction vector 
field is statistically stationary in the long-time limit. This is an apparently paradoxi- 
cal conclusion, in that under the assumption that a Lyapunov exponent is positive, we 
showed that there is not increasing sensitivity to the initial condition in the long-time 
limit. 

For non-steady flows, in the long-time limit the pseudomonodromy matrix becomes 
hyperbolic almost everywhere, and the particles align with its dominant eigenvector. This 
implies that in the long-time limit, the initial condition is forgotten for random velocity 
fields. However there are cases when the pseudomonodromy matrix does not have large 
eigenvalues, so that there is some memory of the initial orientation. In these cases, it is 
usually physically appropriate to assume that the initial particle directions are randomly 
distributed, and to average over a uniform distribution of the initial angle. In these cases 
the typical orientation of the particles is described by an order parameter fleld, C,{r,t). 
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The order parameter is required to characterise the direction field at short times, for any 
type of flow. Also, in the case of recirculating flows, such as those considered here, there 
is no guarant ee that the pseudomonodromy matrix has an eigenvalue which increases 
with time. In iBezug lyy, Wilkinson fcMehlis ( 20091 ) we considered the order parameter 
field for random flows at short times, and showed that this field has true singularities, 
which resemble the core and delta singularities of fingerprints. Their normal forms were 
characterised, and their existence was demonstrated experimentally. 

In this paper we considered the second situation where the order parameter is relevant, 
that of a recirculating flow in the long-time limit. This example turns out to be more 
subtle than the case of random flows at shor t times. This is primaril y beca use it is 
a complement to the theorem proved in .Wilkinson. Bezuglw fc Mehligj ( 20081 ). In this 
case the Lyapunov exponent of the flow is zero, and the pattern formed by the order 
parameter field can depend increasingly sensitively on position as i ^ oo, (witnessed 
by the increasingly tight spirals shown in figure [1]). In order to fully understand the 
evolution of the order parameter, in such cases it is necessary to understand how to 
compute the order parameter of the locally averaged orientation in the long-time limit. 
It is this calculation which is the central achievement of the present work. The result is 
contained in equations (j4.3p and (jS.ip . which show how the matrix defining the inertia 
tensor of the direction distribution is related to the normal form decomposition of the 
transfer matrix Mq. 

Our result on the locally averaged order parameter shows that the expression for the 
quadratic form of the direction inertia tensor (equation (|5.1|) ) has the same structure 
as for the un- averaged case (equation (j3.9p ). This implies that the singularities of the 
averaged order parameter have the same structure as for the un-averaged case. We also 
considered the Poincare indices of the hyperbolic bands, where the particles become 
perfectly aligned. We showed that upon crossing an elliptic band with the topology of 
an annulus, the change of the Poincare index is ±i if the trace of the transfer matrix 
changes sign, and if the sign of tr(Mo) is unchanged. 

Finally we note that our results for recirculating fiows depend upon the aspect ratio of 
the particles (via the parameters ai and a2 = 1 — ai in (|2.2p ). and that in most practical 
applications the particles may not all have the same aspect ratio. This means that the 
boundaries between the elliptic and hyperbolic bands become blurred. In the ideal case, 
these boundaries are only marked by a discontinuity of the order parameter, so that in 
practical applications the boundaries may be very hard to determine. However the zeros 
of the order parameter are much more robust, and their normal forms have the same 
structure even if the rheoscopic suspension has particles which have a disperse aspect 
ratio. 
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